;;; ARRAYS FOR LAYING OUT POISSON DISTRIBUTIONS

(ARRAY P1L1+ FLONUM 80.) 

(ARRAY P0L1- FLONUM 80.) 

(ARRAY P0L0+ FLONUM 80.) 

(ARRAY P1L01- FLONUM 80.) 

(ARRAY P0L00- FLONUM 80.) 

;;; VARIABLES DEFINING LENGTHS OF THE DISTRIBUTIONS
;;; EPSILON IS THE INHIBITION FACTOR

(DECLARE (SPECIAL NP1L1+ NP0L1- NP0L0+ NP1L01- NP0L00- EPSILON)
	 (FIXNUM NP1L1+ NP1L1- NP0L0+ NP1L01- NP0L00-)
	 (FIXNUM L+ L-)
	 (FLONUM F EPSILON)) 

(DECLARE (FIXNUM (ENTER-POISSON-DISTRIBUTION NOTYPE FLONUM) I)
	 (FLONUM VALUE X)) 

;;; SET UP THE DISTRIBUTIONS

(DEFUN ENTER-POISSON-DISTRIBUTION (ARRAYPOINTER MU) 
       (PROG (VALUE X I) 
	     (SETQ VALUE (EXP (-$ MU)) I 0. X 0.0)
	LOOP (STORE (ARRAYCALL FLONUM ARRAYPOINTER I) VALUE)
	     (COND ((AND (> X MU) (< VALUE 1.0E-5)) (RETURN I))
		   (T (SETQ I (1+ I) X (1+$ X))
		      (SETQ VALUE (*$ VALUE (//$ MU X)))
		      (GO LOOP)))))

;;; can do same computations using the more expensive but more
;;; accurate binomial distribution

(DECLARE (FIXNUM N R (ENTER-BINOMIAL-DISTRIBUTION NOTYPE FLONUM FIXNUM))
         (FLONUM P VALUE FACTOR))

(DEFUN ENTER-BINOMIAL-DISTRIBUTION (ARRAYPOINTER P N) 
       (PROG (VALUE FACTOR) 
	     (STORE (ARRAYCALL FLONUM ARRAYPOINTER N) (EXPT P N))
	     (SETQ VALUE (EXPT (-$ 1.0 P) N))
	     (SETQ FACTOR (//$ P (-$ 1.0 P)))
             (STORE (ARRAYCALL FLONUM ARRAYPOINTER 0.) VALUE)
             (DO I 1 (1+ I) (= I N)
                 (SETQ VALUE (*$ VALUE FACTOR (//$ (FLOAT (- N I -1)) (FLOAT I))))
                 (STORE (ARRAYCALL FLONUM ARRAYPOINTER I) VALUE))
             (RETURN (1+ N))))

;;; given L1+, L1- and F, P1, P0 and THETA, derive new P1* and P0*

(DECLARE (SPECIAL L1+ L1- L0+ L01- L00- F P1 P0)
	 (FIXNUM L1+ L1- THETA)
	 (FLONUM F P1 P0 (P1* FIXNUM) (P0* FIXNUM))) 

(SETQ L1+ 13. L1- 12. F (//$ 1.0 7.0) EPSILON 1.0)

(P '(SET L1+ L1- F AND EPSILON DEFAULTS ARE) L1+ L1- F 'AND EPSILON) 

(DEFUN SET-UP-POISSON-DISTRIBUTIONS (PROB1 PROB0) 
       (PROG NIL 
	     ;;; set current probabilities
	     (SETQ P1 PROB1 P0 PROB0)
	     ;;; assumptions about the operator allow these assignments
	     (SETQ L0+ L1+ 
		   L01- (FIX (+$ 0.5 (*$ F (FLOAT L1-)))) 
		   L00- (- L1- L01-))
	     (SETQ NP1L1+
		   (ENTER-POISSON-DISTRIBUTION (GET 'P1L1+
						    'ARRAY)
					       (*$ PROB1
						   (FLOAT L1+))))
	     (SETQ NP0L1-
		   (ENTER-POISSON-DISTRIBUTION (GET 'P0L1-
						    'ARRAY)
					       (*$ PROB0
						   (FLOAT L1-))))
	     (SETQ NP0L0+
		   (ENTER-POISSON-DISTRIBUTION (GET 'P0L0+
						    'ARRAY)
					       (*$ PROB0
						   (FLOAT L0+))))
	     (SETQ NP1L01-
		   (ENTER-POISSON-DISTRIBUTION (GET 'P1L01-
						    'ARRAY)
					       (*$ PROB1
						   (FLOAT L01-))))
	     (SETQ NP0L00-
		   (ENTER-POISSON-DISTRIBUTION (GET 'P0L00-
						    'ARRAY)
					       (*$ PROB0
						   (FLOAT L00-))))
	     (RETURN '(POISSON DISTRIBUTIONS SET UP))))

(DEFUN SET-UP-BINOMIAL-DISTRIBUTIONS (PROB1 PROB0) 
       (PROG NIL 
	     ;;; set current probabilities
	     (SETQ P1 PROB1 P0 PROB0)
	     ;;; assumptions about the operator allow these assignments
	     (SETQ L0+ L1+ 
		   L01- (FIX (+$ 0.5 (*$ F (FLOAT L1-)))) 
		   L00- (- L1- L01-))
	     (SETQ NP1L1+
		   (ENTER-BINOMIAL-DISTRIBUTION (GET 'P1L1+
						    'ARRAY)
					        PROB1
						L1+))
	     (SETQ NP0L1-
		   (ENTER-BINOMIAL-DISTRIBUTION (GET 'P0L1-
						    'ARRAY)
					        PROB0
						L1-))
	     (SETQ NP0L0+
		   (ENTER-BINOMIAL-DISTRIBUTION (GET 'P0L0+
						    'ARRAY)
					        PROB0
						L0+))
	     (SETQ NP1L01-
		   (ENTER-BINOMIAL-DISTRIBUTION (GET 'P1L01-
						    'ARRAY)
					        PROB1
						L01-))
	     (SETQ NP0L00-
		   (ENTER-BINOMIAL-DISTRIBUTION (GET 'P0L00-
						    'ARRAY)
					        PROB0
						L00-))
	     (RETURN '(BINOMIAL DISTRIBUTIONS SET UP)))) 

(DECLARE (FLONUM (PROB-OF-X2 FIXNUM FIXNUM FIXNUM) VALUE THISTERM)
	 (FIXNUM THETA X1 X2 X2LIMIT M1 M2 M3)) 

(DEFUN P1* (THETA) 
       (PROG (M1 M2 X2LIMIT VALUE THISTERM) 
	     ;;; pick up M1 and M2 from the specials
	     (SETQ M1 (1- NP1L1+) M2 (1- NP0L1-))
	     (SETQ VALUE 0.0)
	     (DO X1
		 THETA
		 (1+ X1)
		 (> X1 M1)
		 (SETQ THISTERM 0.0 X2LIMIT (MIN (FIX (//$ (FLOAT (- X1 THETA))
                                                           EPSILON))
                                                  M2))
		 (DO X2
		     0.
		     (1+ X2)
		     (> X2 X2LIMIT)
		     (SETQ THISTERM (+$ THISTERM (P0L1- X2))))
		 (SETQ VALUE (+$ VALUE (*$ (P1L1+ X1) THISTERM))))
	     (RETURN VALUE))) 

(DEFUN P0* (THETA) 
       (PROG (M1 X2LIMIT M2 M3 VALUE THISTERM) 
	     ;;; pick up M1 from the specials
	     (SETQ M1 (1- NP0L0+) M2 (1- NP1L01-) M3 (1- NP0L00-))
	     (SETQ VALUE 0.0)
	     (DO X1
		 THETA
		 (1+ X1)
		 (> X1 M1)
		 (SETQ THISTERM 0.0 X2LIMIT (FIX (//$ (FLOAT (- X1 THETA))
                                                      EPSILON)))
		 (DO X2
		     0.
		     (1+ X2)
		     (> X2 X2LIMIT)
		     (SETQ THISTERM (+$ THISTERM
					(PROB-OF-X2 X2 M2 M3))))
		 (SETQ VALUE (+$ VALUE (*$ (P1L1+ X1) THISTERM))))
	     (RETURN VALUE))) 

(DECLARE (FIXNUM X2 M2 M3 Y1L Y1U Y1) (FLONUM VALUE)) 

(DEFUN PROB-OF-X2 (X2 M2 M3) 
       (PROG (VALUE Y1L Y1U) 
	     (SETQ VALUE 0.0)
	     (SETQ Y1L (MAX 0. (- X2 M3)) Y1U (MIN (1- X2) M2))
	     (DO Y1
		 Y1L
		 (1+ Y1)
		 (> Y1 Y1U)
		 (SETQ VALUE (+$ VALUE
				 (*$ (P1L01- Y1) (P0L00- (- X2 Y1))))))
	     (RETURN VALUE))) 

;;; POISSON DISTRIBUTION

(DECLARE (FIXNUM N (FACT FIXNUM))) 

(DEFUN FACT (N) (COND ((= N 0.) 1.) (T (* N (FACT (1- N)))))) 

(DECLARE (FLONUM (POISSON FIXNUM FLONUM))) 

(DEFUN POISSON (N MU) 
       (COND ((= N 0.) (EXP (-$ MU)))
	     (T (*$ (POISSON (1- N) MU) (//$ MU (FLOAT N)))))) 

(DECLARE (FLONUM A B ANS (BINOM FIXNUM FIXNUM))) 

(DEFUN BINOM (N R) 
       (PROG (B ANS) 
	     (SETQ ANS 1.0)
	     (SETQ B (FLOAT N))
	     (DO A
		 (FLOAT R)
		 (1-$ A)
		 (= A 0.0)
		 (SETQ ANS (//$ (*$ ANS B) A))
		 (SETQ B (1-$ B)))
	     (RETURN ANS))) 

(DECLARE (FLONUM X (BINOMIAL FLONUM FIXNUM FIXNUM)) (FIXNUM N R)) 

(DEFUN BINOMIAL (X R N) 
       (*$ (BINOM N R) (EXPT X R) (EXPT (-$ 1.0 X) (- N R)))) 
βββ